#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

sig.effects(measure = "FA", atlas = "glasser", dataset = "hbn")

smoothcentered_FA_glasser_hbn.plot <- merge((smoothcentered_FA_glasser_hbn %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hbn %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_hbn.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.017, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(6, 8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Developmental_Trajectories_HBN.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.6)

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hbn, by = c("tract", "label", "orig_parcelname"), suffixes = c("_pnc", "_hbn"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hbn)

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hbn, perm.id = perm.id.full, corr.type = "pearson") 

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hbn)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.051, 0.15)) +
  labs(x="\nPNC", y="HBN\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/FA_Ageeffect_correlation_PNCHBN.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.6)

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hbn)

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hbn, perm.id = perm.id.full, corr.type = "pearson") 

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

gameffects_neurosynth_hbn <- merge(gameffects_FA_glasser_hbn, neurosynth.terms, by = "label")

devpattern.neurosynth.hbn <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hbn, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))

devpattern.neurosynth.hbn <- rbind(slice_max(devpattern.neurosynth.hbn, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hbn, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.hbn

ggplot(devpattern.neurosynth.hbn, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 1.8, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.6,0.62), breaks = c(-0.6, -0.3, 0, 0.3, 0.6)) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_blank(),
   panel.grid.minor.y = element_blank(),
   panel.grid.minor.x = element_blank(),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Neurosynth_Maturationdecoding_HBN.pdf", device = "pdf", dpi = 500, width = 1.9, height = 3.1)

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))

devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hbn, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)

derivatives_FA_glasser_hbn$significant.derivative[derivatives_FA_glasser_hbn$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hbn, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(linewidth = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

cor.test(gameffects_FA_glasser_hbn$SA.axis, gameffects_FA_glasser_hbn$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hbn, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(gameffects_FA_glasser_hbn, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(11, 20), breaks = c(12, 14, 16, 18, 20)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/SAaxis_agematuration_hbn.pdf", device = "pdf", dpi = 500, width = 1.9, height = 1.6)

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)

#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)

EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 

EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)

df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

plasticity.dev.hbn <- left_join(plasticity.dev, gameffects_FA_glasser_hbn, by = c("label", "orig_parcelname", "tract"))

cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$EI.ageslope, method = c("spearman"), exact = F)

EI.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hbn, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#fcd16d") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$myelin.annualroc, method = c("spearman"), exact = F)

myelin.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hbn, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#F1A45F") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$amplitude.agedecline, method = c("spearman"), exact = F)

amplitude.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

ggplot(plasticity.dev.hbn, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#7B2C80") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

plasticity.ps <- c(EI.thalamus.hbn, myelin.thalamus.hbn, amplitude.thalamus.hbn)
plasticity.ps

p.adjust(plasticity.ps, method = c("fdr"))

plasticity.results <- data.frame(map = factor(), corr = double())
plasticity.results <- plasticity.results %>% add_row(map = "E/I Ratio", corr = 0.5698692)
plasticity.results <- plasticity.results %>% add_row(map = "T1/T2 Ratio", corr = -0.5783632)
plasticity.results <- plasticity.results %>% add_row(map = "BOLD Amplitude", corr = 0.6852288)
plasticity.results$map <- factor(plasticity.results$map, ordered = TRUE, levels = c( "E/I Ratio", "T1/T2 Ratio", "BOLD Amplitude"))

ggplot(plasticity.results, aes(x = map, y = corr)) + 
  geom_segment(aes(x = map, xend = map,  yend =  
                     corr, y = 0), linewidth = 0.2, color = "#989898") +
  geom_point(size = 2, alpha = 0.85, color = "#FCAB6A", fill = "#FCAB6A", shape = 22) +
  labs(x="Cortical Development Map") +
  labs(y="\nThalamocortical Correlation (r)") +
  theme_classic() +
  scale_y_continuous(breaks = c(-0.6, -0.3, 0, 0.3, 0.6), limits = c(-0.62, 0.7)) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.text.x = element_text(vjust = 0.7, angle = 20),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Plasticitymap_correlation_HBN.pdf", device = "pdf", dpi = 500, width = 1.9, height = 1.58)

generalSES_maineffects_FA_glasser_hbn$significant <- p.adjust(generalSES_maineffects_FA_glasser_hbn$GAM.cov.pvalue, method = c("fdr")) < 0.05
sum(generalSES_maineffects_FA_glasser_hbn$significant)

sum(generalSES_maineffects_FA_glasser_hbn$significant)/length(generalSES_maineffects_FA_glasser_hbn$significant)

sigeffects.df <- generalSES_maineffects_FA_glasser_hbn %>% filter(significant == TRUE)

(sum(sigeffects.df$significant > 0)/nrow(sigeffects.df))*100

spin.data <- left_join(spin.parcels, sigeffects.df, by = "tract") #merged data for spin tests, with all 360 parcels in the expected right hemisphere --> left hemisphere order (matching glasser.coords)
spin.data <- merge(spin.data, S.A.axis, by = c("label", "orig_parcelname", "tract"), sort =F)

cor.test(spin.data$GAM.cov.tvalue, spin.data$SA.axis, method = c("spearman"))

perm.sphere.p(x = spin.data$SA.axis, y = spin.data$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") 

ggplot() +
  geom_point(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue, color = SA.axis), size = 0.5) +
  geom_smooth(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue), method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  labs(x="\nS-A axis", y="Environment effect (t)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/EnvEffect_SArank_correlationplot_HBN.pdf", device = "pdf", width = 1.9, height = 1.6)

#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>% 
                             group_by(SA.axis.bin) %>% 
                             do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
                             unnest(cols = c("mean.tvalue.significant"))

envSES_effects_SAaxisbins

#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)

## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix 
x <- spin.data$SA.axis.bin #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full

## spin the empirical data 
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
  for (r in 1:nperm) {
    for (i in 1:nroi) {
      x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
      y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
    }
  }

## compute mean t-value for each S-A axis bin using spatially permuted x inputs 
### set up spun (x) and empirical (y) df
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm.empiricaly <- cbind(y, x.spatialperm) %>% as.data.frame()
colnames(x.spatialperm.empiricaly)[1] <- c("empirical.y")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedx <- sapply(x.spatialperm.empiricaly[2:ncol(x.spatialperm.empiricaly)], #for all permuted bins
                                              function(p){
                                                x.spatialperm.empiricaly %>% group_by(as.character(p)) %>% #group by permuted S-A axis bins
                                                do(mean.tvalue.significant = mean(.$empirical.y, na.rm = TRUE)) %>% #calculate mean t-value
                                                unnest(cols = c("mean.tvalue.significant")) %>% t()}) %>% 
                                                t() %>% as.data.frame()
                                                 
envSES_effects_SAaxisbins_permutedx <- envSES_effects_SAaxisbins_permutedx %>% mutate_if(is.character, as.numeric) #make it numeric
envSES_effects_SAaxisbins_permutedx <-  envSES_effects_SAaxisbins_permutedx[ , !c(TRUE,FALSE) ] 
colnames(envSES_effects_SAaxisbins_permutedx) <- c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue") 

## compute mean t-value for each S-A axis bin using spatially permuted y inputs 
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
                                       dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
                                       select(-empirical.x) %>% t() %>% as.data.frame() %>%
                                       set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue"))

histogram.data <- rbind(envSES_effects_SAaxisbins_permutedx, envSES_effects_SAaxisbins_permutedy)
saveRDS(histogram.data, "/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.HBN.rds")

histogram.data <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.HBN.rds")
envSES_effects_SAaxisbins_permutedx <- histogram.data[1:10000,]
envSES_effects_SAaxisbins_permutedy <- histogram.data[10001:20000,]
nperm = dim(perm.id.full)[2]

#Function to plot the null distribution of environment effect t-values for a given S-A axis bin and test whether the true t-value is significantly smaller or greater in magnitude than the null (permutation based p)
enveffects_SAbin_enrichment <- function(quintile, enrichment_type){
  
  SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.bin == quintile) %>% select(mean.tvalue.significant)
  
  #Histogram
  SA.colorvector <- c("#FFC228", "#F9D58C", "#e9dcf2", "#BE93C4", "#6F1382")
  
  null.tvalue.histogram <- histogram.data %>% select(sprintf("SA.bin%s.tvalue", quintile)) %>%
    ggplot(aes(x = .[,1])) + geom_histogram(fill = c(SA.colorvector)[quintile]) +
    theme_classic() +
    geom_vline(xintercept =  SAbin.tvalue.empirical$mean.tvalue.significant, linewidth = 0.25) +
    xlab("Environment effect (t-value)") +
    ylab("Number of nulls") +
    scale_x_continuous(breaks = c(2, 3, 4), limits = c(2, 4)) +
    theme(
     legend.position = "none", 
     axis.text.x = element_text(size = 5, family = "Arial", color = c("black")),
     axis.text.y = element_blank(),
     axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
     axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
     axis.line.x = element_line(linewidth = .2), 
     axis.line.y = element_blank(),
     axis.ticks = element_line(linewidth = .2))
    
  #Permutation-based p
    if(enrichment_type == "smaller"){
      p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- (p.perm.xy+p.perm.yx)/2
    }
  
  if(enrichment_type == "larger"){
      p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- p.perm.xy
      p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- (p.perm.xy+p.perm.yx)/2
  }
  
  permutation.results <- list(null.tvalue.histogram, perm.p)
  return(permutation.results)
}

enveffects_SAbin_enrichment(quintile = 1, enrichment_type = "smaller")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile1_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

enveffects_SAbin_enrichment(quintile = 2, enrichment_type = "smaller")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile2_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

enveffects_SAbin_enrichment(quintile = 3, enrichment_type = "smaller")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile3_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

enveffects_SAbin_enrichment(quintile = 4, enrichment_type = "larger")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile4_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

enveffects_SAbin_enrichment(quintile = 5, enrichment_type = "larger")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile5_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)